
const char *short_source =
"__kernel void Triad(__global const float *memA, __global const float *memB, __global float *memC)""\n"
"{""\n"
"    int gid = get_global_id(0);""\n"
"    memC[gid] = memA[gid] + memB[gid];""\n"
"}""\n";


const char *long_source =
"__kernel void""\n"
"uniformAdd(__global float *g_data,""\n"
"           __global float *uniforms,""\n"
"           int n,""\n"
"           int blockOffset,""\n"
"           int baseIndex)""\n"
"{""\n"
"    float uni = 0.0f;""\n"
"\n"
"    uni = uniforms[get_group_id(0) + blockOffset];""\n"
"    unsigned int address = (get_group_id(0) * (get_local_size(0) << 1)) +""\n"
"                           baseIndex + get_local_id(0);""\n"
"\n"
"    barrier(CLK_LOCAL_MEM_FENCE);""\n"
"\n"
"    g_data[address] += uni;""\n"
"    if (get_local_id(0) + get_local_size(0) < n)""\n"
"    {""\n"
"        g_data[address + get_local_size(0)] +=  uni;""\n"
"    }""\n"
"}""\n"
"\n"
"__kernel void""\n"
"scan(__global float *g_odata,""\n"
"     __global float *g_idata,""\n"
"     __global float *g_blockSums,""\n"
"     int n,""\n"
"     int blockIndex,""\n"
"     int baseIndex,""\n"
"     int storeSum,""\n"
"     __local float *s_data)""\n"
"{""\n"
"    int ai, bi;""\n"
"    int mem_ai, mem_bi;""\n"
"    int bIndex;""\n"
"\n"
"    // load data into shared memory""\n"
"    if (baseIndex == 0)""\n"
"    {""\n"
"        bIndex = get_group_id(0) * (get_local_size(0) << 1);""\n"
"    }""\n"
"    else""\n"
"    {""\n"
"        bIndex = baseIndex;""\n"
"    }""\n"
"\n"
"    int thid = get_local_id(0);""\n"
"    mem_ai = bIndex + thid;""\n"
"    mem_bi = mem_ai + get_local_size(0);""\n"
"\n"
"    ai = thid;""\n"
"    bi = thid + get_local_size(0);""\n"
"\n"
"    // Cache the computational window in shared memory""\n"
"    // pad values beyond n with zeros""\n"
"    s_data[ai] = g_idata[mem_ai];""\n"
"    if (bi < n)""\n"
"    {""\n"
"        s_data[bi] = g_idata[mem_bi];""\n"
"    }""\n"
"    else""\n"
"    {""\n"
"        s_data[bi] = 0.0f;""\n"
"    }""\n"
"\n"
"    unsigned int stride = 1;""\n"
"\n"
"    // build the sum in place up the tree""\n"
"    for (int d = get_local_size(0); d > 0; d >>= 1)""\n"
"    {""\n"
"        barrier(CLK_LOCAL_MEM_FENCE);""\n"
"        if (thid < d)""\n"
"        {""\n"
"            int i  = 2 * stride * thid;""\n"
"            int aii = i + stride - 1;""\n"
"            int bii = aii + stride;""\n"
"\n"
"            s_data[bii] += s_data[aii];""\n"
"        }""\n"
"        stride *= 2;""\n"
"    }""\n"
"\n"
"    bIndex = (blockIndex == 0) ? get_group_id(0) : blockIndex;""\n"
"\n"
"    if (get_local_id(0) == 0)""\n"
"    {""\n"
"        int index = (get_local_size(0) << 1) - 1;""\n"
"\n"
"        if (storeSum == 1)""\n"
"        {""\n"
"            // write this block's total sum to the corresponding""\n"
"            // index in the blockSums array""\n"
"            g_blockSums[bIndex] = s_data[index];""\n"
"        }""\n"
"\n"
"        // zero the last element in the scan so it will propagate""\n"
"        // back to the front""\n"
"        s_data[index] = 0;""\n"
"    }""\n"
"\n"
"    barrier(CLK_LOCAL_MEM_FENCE);""\n"
"\n"
"    // traverse down the tree building the scan in place""\n"
"    for (int d = 1; d <= get_local_size(0); d *= 2)""\n"
"    {""\n"
"        stride >>= 1;""\n"
"        barrier(CLK_LOCAL_MEM_FENCE);""\n"
"\n"
"        if (thid < d)""\n"
"        {""\n"
"            int i  = 2 * stride * thid;""\n"
"            int aii = i + stride - 1;""\n"
"            int bii = aii + stride;""\n"
"\n"
"            float t  = s_data[aii];""\n"
"            s_data[aii] = s_data[bii];""\n"
"            s_data[bii] += t;""\n"
"        }""\n"
"    }""\n"
"\n"
"    barrier(CLK_LOCAL_MEM_FENCE);""\n"
"\n"
"    // write results to global memory""\n"
"    g_odata[mem_ai] = s_data[ai];""\n"
"    if (bi < n)""\n"
"    {""\n"
"        g_odata[mem_bi] = s_data[bi];""\n"
"    }""\n"
"}""\n";

const char * const matrix_kernels = 
"// file automatically generated - do not edit!\n"
"// matrix-matrix multiplication C = A^T * B^T\n"
"// matrix layouts: C...col_major, A...row_major, B...col_major\n"
"__kernel void prod_TT(\n"
"          __global const float * A,\n"
"          unsigned int A_rows,\n"
"          unsigned int A_cols,\n"
"          unsigned int A_internal_rows,\n"
"          unsigned int A_internal_cols,\n"
"          __global const float * B,  \n"
"          unsigned int B_rows,\n"
"          unsigned int B_cols,\n"
"          unsigned int B_internal_rows,\n"
"          unsigned int B_internal_cols,\n"
"          __global float * C,\n"
"          unsigned int C_rows,\n"
"          unsigned int C_cols,\n"
"          unsigned int C_internal_rows,\n"
"          unsigned int C_internal_cols,\n"
"          __local float * bufA,\n"
"          __local float * bufB) \n"
"{ \n"
"  int block_size = get_local_size(0);\n"
"  int row_block_id = get_group_id(0);\n"
"  int col_block_id = get_group_id(1);\n"
"  int row_thread_id = get_local_id(0);\n"
"  int col_thread_id = get_local_id(1);\n"
"  int aBegin = row_block_id * block_size;\n"
"  int aStep = block_size * A_internal_cols;\n"
"  int bBegin = col_block_id * block_size;\n"
"  int bStep = block_size * B_internal_rows;\n"
"  int block_num = A_rows / block_size;\n"
"  if (block_num * block_size != A_rows)\n"
"    ++block_num;\n"
"  float Csub = 0;\n"
"  int aOffset = row_thread_id + col_thread_id * A_internal_cols;\n"
"  int bOffset = row_thread_id * B_internal_rows + col_thread_id;\n"
"  for (int block = 0;\n"
"           block < block_num;\n"
"           ++block)\n"
"  {\n"
"    if (block * block_size + col_thread_id < A_rows && get_global_id(0) < A_cols)\n"
"      bufA[row_thread_id * block_size + col_thread_id] = A[aBegin + aOffset];\n"
"    else\n"
"      bufA[row_thread_id * block_size + col_thread_id] = 0;\n"
"    if ( (block * block_size + row_thread_id < B_cols) && get_global_id(1) < B_rows )\n"
"      bufB[row_thread_id * block_size + col_thread_id] = B[bBegin + bOffset]; \n"
"    else\n"
"      bufB[row_thread_id * block_size + col_thread_id] = 0;\n"
"    barrier(CLK_LOCAL_MEM_FENCE);\n"
"    for (int k = 0; k < block_size; ++k)\n"
"      Csub += bufA[row_thread_id * block_size + k] * bufB[k * block_size + col_thread_id];\n"
"    barrier(CLK_LOCAL_MEM_FENCE);\n"
"    aBegin += aStep;\n"
"    bBegin += bStep;\n"
"  }\n"
"  if (get_global_id(0) < A_cols && get_global_id(1) < B_rows)\n"
"    C[get_global_id(0) + get_global_id(1) * C_internal_rows] = Csub;\n"
"}\n"
"// file automatically generated - do not edit!\n"
"// matrix-matrix multiplication C = A^T * B\n"
"// matrix layouts: C...col_major, A...row_major, B...col_major\n"
"__kernel void prod_TA(\n"
"          __global const float * A,\n"
"          unsigned int A_rows,\n"
"          unsigned int A_cols,\n"
"          unsigned int A_internal_rows,\n"
"          unsigned int A_internal_cols,\n"
"          __global const float * B,  \n"
"          unsigned int B_rows,\n"
"          unsigned int B_cols,\n"
"          unsigned int B_internal_rows,\n"
"          unsigned int B_internal_cols,\n"
"          __global float * C,\n"
"          unsigned int C_rows,\n"
"          unsigned int C_cols,\n"
"          unsigned int C_internal_rows,\n"
"          unsigned int C_internal_cols,\n"
"          __local float * bufA,\n"
"          __local float * bufB) \n"
"{ \n"
"  int block_size = get_local_size(0);\n"
"  int row_block_id = get_group_id(0);\n"
"  int col_block_id = get_group_id(1);\n"
"  int row_thread_id = get_local_id(0);\n"
"  int col_thread_id = get_local_id(1);\n"
"  int aBegin = row_block_id * block_size;\n"
"  int aStep = block_size * A_internal_cols;\n"
"  int bBegin = col_block_id * block_size * B_internal_rows;\n"
"  int bStep = block_size;\n"
"  int block_num = A_rows / block_size;\n"
"  if (block_num * block_size != A_rows)\n"
"    ++block_num;\n"
"  float Csub = 0;\n"
"  int aOffset = row_thread_id + col_thread_id * A_internal_cols;\n"
"  int bOffset = row_thread_id + col_thread_id * B_internal_rows;\n"
"  for (int block = 0;\n"
"           block < block_num;\n"
"           ++block)\n"
"  {\n"
"    if (block * block_size + col_thread_id < A_rows && get_global_id(0) < A_cols)\n"
"      bufA[row_thread_id * block_size + col_thread_id] = A[aBegin + aOffset];\n"
"    else\n"
"      bufA[row_thread_id * block_size + col_thread_id] = 0;\n"
"    if ( (block * block_size + row_thread_id < B_rows) && get_global_id(1) < B_cols )\n"
"      bufB[row_thread_id * block_size + col_thread_id] = B[bBegin + bOffset]; \n"
"    else\n"
"      bufB[row_thread_id * block_size + col_thread_id] = 0;\n"
"    barrier(CLK_LOCAL_MEM_FENCE);\n"
"    for (int k = 0; k < block_size; ++k)\n"
"      Csub += bufA[row_thread_id * block_size + k] * bufB[k * block_size + col_thread_id];\n"
"    barrier(CLK_LOCAL_MEM_FENCE);\n"
"    aBegin += aStep;\n"
"    bBegin += bStep;\n"
"  }\n"
"  if (get_global_id(0) < A_cols && get_global_id(1) < B_cols)\n"
"    C[get_global_id(0) + get_global_id(1) * C_internal_rows] = Csub;\n"
"}\n"
"// file automatically generated - do not edit!\n"
"// matrix-matrix multiplication C = A * B\n"
"// matrix layouts: C...col_major, A...row_major, B...col_major\n"
"__kernel void prod_AA(\n"
"          __global const float * A,\n"
"          unsigned int A_rows,\n"
"          unsigned int A_cols,\n"
"          unsigned int A_internal_rows,\n"
"          unsigned int A_internal_cols,\n"
"          __global const float * B,  \n"
"          unsigned int B_rows,\n"
"          unsigned int B_cols,\n"
"          unsigned int B_internal_rows,\n"
"          unsigned int B_internal_cols,\n"
"          __global float * C,\n"
"          unsigned int C_rows,\n"
"          unsigned int C_cols,\n"
"          unsigned int C_internal_rows,\n"
"          unsigned int C_internal_cols,\n"
"          __local float * bufA,\n"
"          __local float * bufB) \n"
"{ \n"
"  int block_size = get_local_size(0);\n"
"  int row_block_id = get_group_id(0);\n"
"  int col_block_id = get_group_id(1);\n"
"  int row_thread_id = get_local_id(0);\n"
"  int col_thread_id = get_local_id(1);\n"
"  int aBegin = row_block_id * block_size * A_internal_cols;\n"
"  int aStep = block_size;\n"
"  int bBegin = col_block_id * block_size * B_internal_rows;\n"
"  int bStep = block_size;\n"
"  int block_num = A_cols / block_size;\n"
"  if (block_num * block_size != A_cols)\n"
"    ++block_num;\n"
"  float Csub = 0;\n"
"  int aOffset = row_thread_id * A_internal_cols + col_thread_id;\n"
"  int bOffset = row_thread_id + col_thread_id * B_internal_rows;\n"
"  for (int block = 0;\n"
"           block < block_num;\n"
"           ++block)\n"
"  {\n"
"    if (block * block_size + col_thread_id < A_cols && get_global_id(0) < A_rows)\n"
"      bufA[row_thread_id * block_size + col_thread_id] = A[aBegin + aOffset];\n"
"    else\n"
"      bufA[row_thread_id * block_size + col_thread_id] = 0;\n"
"    if ( (block * block_size + row_thread_id < B_rows) && get_global_id(1) < B_cols )\n"
"      bufB[row_thread_id * block_size + col_thread_id] = B[bBegin + bOffset]; \n"
"    else\n"
"      bufB[row_thread_id * block_size + col_thread_id] = 0;\n"
"    barrier(CLK_LOCAL_MEM_FENCE);\n"
"    for (int k = 0; k < block_size; ++k)\n"
"      Csub += bufA[row_thread_id * block_size + k] * bufB[k * block_size + col_thread_id];\n"
"    barrier(CLK_LOCAL_MEM_FENCE);\n"
"    aBegin += aStep;\n"
"    bBegin += bStep;\n"
"  }\n"
"  if (get_global_id(0) < A_rows && get_global_id(1) < B_cols)\n"
"    C[get_global_id(0) + get_global_id(1) * C_internal_rows] = Csub;\n"
"}\n"
"// file automatically generated - do not edit!\n"
"// matrix-matrix multiplication C = A * B^T\n"
"// matrix layouts: C...col_major, A...row_major, B...col_major\n"
"__kernel void prod_AT(\n"
"          __global const float * A,\n"
"          unsigned int A_rows,\n"
"          unsigned int A_cols,\n"
"          unsigned int A_internal_rows,\n"
"          unsigned int A_internal_cols,\n"
"          __global const float * B,  \n"
"          unsigned int B_rows,\n"
"          unsigned int B_cols,\n"
"          unsigned int B_internal_rows,\n"
"          unsigned int B_internal_cols,\n"
"          __global float * C,\n"
"          unsigned int C_rows,\n"
"          unsigned int C_cols,\n"
"          unsigned int C_internal_rows,\n"
"          unsigned int C_internal_cols,\n"
"          __local float * bufA,\n"
"          __local float * bufB) \n"
"{ \n"
"  int block_size = get_local_size(0);\n"
"  int row_block_id = get_group_id(0);\n"
"  int col_block_id = get_group_id(1);\n"
"  int row_thread_id = get_local_id(0);\n"
"  int col_thread_id = get_local_id(1);\n"
"  int aBegin = row_block_id * block_size * A_internal_cols;\n"
"  int aStep = block_size;\n"
"  int bBegin = col_block_id * block_size;\n"
"  int bStep = block_size * B_internal_rows;\n"
"  int block_num = A_cols / block_size;\n"
"  if (block_num * block_size != A_cols)\n"
"    ++block_num;\n"
"  float Csub = 0;\n"
"  int aOffset = row_thread_id * A_internal_cols + col_thread_id;\n"
"  int bOffset = row_thread_id * B_internal_rows + col_thread_id;\n"
"  for (int block = 0;\n"
"           block < block_num;\n"
"           ++block)\n"
"  {\n"
"    if (block * block_size + col_thread_id < A_cols && get_global_id(0) < A_rows)\n"
"      bufA[row_thread_id * block_size + col_thread_id] = A[aBegin + aOffset];\n"
"    else\n"
"      bufA[row_thread_id * block_size + col_thread_id] = 0;\n"
"    if ( (block * block_size + row_thread_id < B_cols) && get_global_id(1) < B_rows )\n"
"      bufB[row_thread_id * block_size + col_thread_id] = B[bBegin + bOffset]; \n"
"    else\n"
"      bufB[row_thread_id * block_size + col_thread_id] = 0;\n"
"    barrier(CLK_LOCAL_MEM_FENCE);\n"
"    for (int k = 0; k < block_size; ++k)\n"
"      Csub += bufA[row_thread_id * block_size + k] * bufB[k * block_size + col_thread_id];\n"
"    barrier(CLK_LOCAL_MEM_FENCE);\n"
"    aBegin += aStep;\n"
"    bBegin += bStep;\n"
"  }\n"
"  if (get_global_id(0) < A_rows && get_global_id(1) < B_rows)\n"
"    C[get_global_id(0) + get_global_id(1) * C_internal_rows] = Csub;\n"
"}\n"
"__kernel void unit_lower_triangular_substitute_inplace(\n"
"          __global const float * matrix,\n"
"          unsigned int matrix_rows,\n"
"          unsigned int matrix_cols,\n"
"          unsigned int matrix_internal_rows,\n"
"          unsigned int matrix_internal_cols,\n"
"          __global float * vector)\n"
"{\n"
"  float temp;\n"
"  for (int row = 0; row < matrix_rows; ++row)\n"
"  {\n"
"    barrier(CLK_GLOBAL_MEM_FENCE);\n"
"    temp = vector[row];\n"
"    for  (int elim = row + get_global_id(0) + 1; elim < matrix_rows; elim += get_global_size(0))\n"
"      vector[elim] -= temp * matrix[row * matrix_internal_rows + elim];\n"
"  }\n"
"}\n"
"__kernel void inplace_sub(\n"
"          __global float * vec1,\n"
"          __global const float * vec2,\n"
"          unsigned int size) \n"
"{ \n"
"  for (unsigned int i = get_global_id(0); i < size; i += get_global_size(0))\n"
"    vec1[i] -= vec2[i];\n"
"}\n"
"__kernel void lower_triangular_substitute_inplace(\n"
"          __global const float * matrix,\n"
"          unsigned int matrix_rows,\n"
"          unsigned int matrix_cols,\n"
"          unsigned int matrix_internal_rows,\n"
"          unsigned int matrix_internal_cols,\n"
"          __global float * vector)\n"
"{\n"
"  float temp;\n"
"  for (int row = 0; row < matrix_rows; ++row)\n"
"  {\n"
"    barrier(CLK_GLOBAL_MEM_FENCE);\n"
"    if (get_global_id(0) == 0)\n"
"      vector[row] /= matrix[row+row*matrix_internal_rows];\n"
"    barrier(CLK_GLOBAL_MEM_FENCE);\n"
"    temp = vector[row];\n"
"    for  (int elim = row + get_global_id(0) + 1; elim < matrix_rows; elim += get_global_size(0))\n"
"      vector[elim] -= temp * matrix[row * matrix_internal_rows + elim];\n"
"  }\n"
"}\n"
"__kernel void trans_vec_mul(\n"
"          __global const float * matrix,\n"
"          unsigned int matrix_rows,\n"
"          unsigned int matrix_cols,\n"
"          unsigned int matrix_internal_rows,\n"
"          unsigned int matrix_internal_cols,\n"
"          __global const float * vector,  \n"
"          __global float * result) \n"
"{ \n"
"  //row and col indicate indices within transposed matrix\n"
"  for (unsigned int row = get_global_id(0); row < matrix_cols; row += get_global_size(0))\n"
"  {\n"
"    float dot_prod2 = 0.0f;\n"
"    for (unsigned int col = 0; col < matrix_rows; ++col)\n"
"      dot_prod2 += matrix[row * matrix_internal_rows + col] * vector[col];\n"
"    result[row] = dot_prod2;\n"
"  }\n"
"}\n"
"//perform a rank-1 update of the matrix, i.e. A += x * x^T\n"
"__kernel void rank1_update(\n"
"          __global float * matrix,\n"
"          unsigned int matrix_rows,\n"
"          unsigned int matrix_cols,\n"
"          unsigned int matrix_internal_rows,\n"
"          unsigned int matrix_internal_cols,\n"
"          __global const float * vector1,  \n"
"          __global const float * vector2) \n"
"{ \n"
"  float tmp;\n"
"  for (unsigned int row= get_global_id(0); row < matrix_rows; row += get_global_size(0))\n"
"  {\n"
"    tmp = vector1[row];\n"
"    for (unsigned int col = 0; col < matrix_cols; ++col)\n"
"      matrix[row + col * matrix_internal_rows] += tmp * vector2[col];\n"
"  }\n"
"}\n"
"__kernel void sub(\n"
"          __global const float * vec1,\n"
"          __global const float * vec2, \n"
"          __global float * result,\n"
"          unsigned int size)\n"
"{ \n"
"  for (unsigned int i = get_global_id(0); i < size; i += get_global_size(0))\n"
"    result[i] = vec1[i] - vec2[i];\n"
"}\n"
"//transposed lower triangular matrix\n"
"__kernel void trans_unit_upper_triangular_substitute_inplace(\n"
"          __global const float * matrix, \n"
"          unsigned int matrix_rows,\n"
"          unsigned int matrix_cols,\n"
"          unsigned int matrix_internal_rows,\n"
"          unsigned int matrix_internal_cols,\n"
"          __global float * vector) \n"
"{ \n"
"  float temp; \n"
"  for (int row = matrix_rows-1; row > -1; --row) \n"
"  { \n"
"    barrier(CLK_GLOBAL_MEM_FENCE); \n"
"    \n"
"    temp = vector[row]; \n"
"    //eliminate column with index 'row' in parallel: \n"
"    for  (int elim = get_global_id(0); elim < row; elim += get_global_size(0)) \n"
"      vector[elim] -= temp * matrix[row + elim  * matrix_internal_rows]; \n"
"  } \n"
"   \n"
"}\n"
"__kernel void lu_factorize(\n"
"          __global float * matrix,\n"
"          unsigned int matrix_rows,\n"
"          unsigned int matrix_cols,\n"
"          unsigned int matrix_internal_rows,\n"
"          unsigned int matrix_internal_cols) \n"
"{ \n"
"  float temp;\n"
"  for (unsigned int i=1; i<matrix_rows; ++i)\n"
"  {\n"
"    for (unsigned int k=0; k<i; ++k)\n"
"    {\n"
"      if (get_global_id(0) == 0)\n"
"        matrix[i + k*matrix_internal_rows] /= matrix[k + k*matrix_internal_rows];\n"
"      barrier(CLK_GLOBAL_MEM_FENCE);\n"
"      temp = matrix[i + k*matrix_internal_rows];\n"
"      \n"
"      //parallel subtraction:\n"
"      for (unsigned int j=k+1 + get_global_id(0); j<matrix_cols; j += get_global_size(0))\n"
"        matrix[i + j*matrix_internal_rows] -= temp * matrix[k + j*matrix_internal_rows];\n"
"    }\n"
"  }\n"
"} \n"
"__kernel void add(\n"
"          __global const float * vec1,\n"
"          __global const float * vec2, \n"
"          __global float * result,\n"
"          unsigned int size) \n"
"{ \n"
"  for (unsigned int i = get_global_id(0); i < size; i += get_global_size(0))\n"
"    result[i] = vec1[i] + vec2[i];\n"
"}\n"
"__kernel void vec_mul(\n"
"          __global const float * matrix,\n"
"          unsigned int matrix_rows,\n"
"          unsigned int matrix_cols,\n"
"          unsigned int matrix_internal_rows,\n"
"          unsigned int matrix_internal_cols,\n"
"          __global const float * vector,  \n"
"          __global float * result) \n"
"{ \n"
"  for (unsigned int row = get_global_id(0); row < matrix_rows; row += get_global_size(0))\n"
"  {\n"
"    float dot_prod = 0.0f;\n"
"    for (unsigned int col = 0; col < matrix_cols; ++col)\n"
"      dot_prod += matrix[row + col*matrix_internal_rows] * vector[col];\n"
"    result[row] = dot_prod;\n"
"  }\n"
"}\n"
"__kernel void trans_lower_triangular_substitute_inplace(\n"
"          __global const float * matrix,\n"
"          unsigned int matrix_rows,\n"
"          unsigned int matrix_cols,\n"
"          unsigned int matrix_internal_rows,\n"
"          unsigned int matrix_internal_cols,\n"
"          __global float * vector)\n"
"{\n"
"  float temp;\n"
"  for (int row = 0; row < matrix_rows; ++row)\n"
"  {\n"
"    barrier(CLK_GLOBAL_MEM_FENCE);\n"
"    if (get_global_id(0) == 0)\n"
"      vector[row] /= matrix[row+row*matrix_internal_rows];\n"
"    barrier(CLK_GLOBAL_MEM_FENCE);\n"
"    temp = vector[row];\n"
"    for  (int elim = row + get_global_id(0) + 1; elim < matrix_rows; elim += get_global_size(0))\n"
"      vector[elim] -= temp * matrix[elim * matrix_internal_rows + row];\n"
"  }\n"
"}\n"
"__kernel void inplace_divide(\n"
"          __global float * vec,\n"
"          __global const float * fac,  //note: CPU variant is mapped to prod_scalar\n"
"          unsigned int size) \n"
"{ \n"
"  float factor = *fac;\n"
"  for (unsigned int i = get_global_id(0); i < size; i += get_global_size(0))\n"
"    vec[i] /= factor;\n"
"}\n"
"//transposed lower triangular matrix\n"
"__kernel void trans_upper_triangular_substitute_inplace(\n"
"          __global const float * matrix, \n"
"          unsigned int matrix_rows,\n"
"          unsigned int matrix_cols,\n"
"          unsigned int matrix_internal_rows,\n"
"          unsigned int matrix_internal_cols,\n"
"          __global float * vector) \n"
"{ \n"
"  float temp; \n"
"  for (int row = matrix_rows-1; row > -1; --row) \n"
"  { \n"
"    barrier(CLK_GLOBAL_MEM_FENCE); \n"
"    if (get_global_id(0) == 0) \n"
"      vector[row] /= matrix[row + row*matrix_internal_rows]; \n"
" \n"
"    barrier(CLK_GLOBAL_MEM_FENCE); \n"
"    temp = vector[row]; \n"
"    //eliminate column with index 'row' in parallel: \n"
"    for  (int elim = get_global_id(0); elim < row; elim += get_global_size(0)) \n"
"      vector[elim] -= temp * matrix[row + elim * matrix_internal_rows]; \n"
"  } \n"
"   \n"
"}\n"
"__kernel void unit_upper_triangular_substitute_inplace( \n"
"          __global const float * matrix, \n"
"          unsigned int matrix_rows,\n"
"          unsigned int matrix_cols,\n"
"          unsigned int matrix_internal_rows,\n"
"          unsigned int matrix_internal_cols,\n"
"          __global float * vector) \n"
"{ \n"
"  float temp; \n"
"  for (int row = matrix_rows-1; row > -1; --row) \n"
"  { \n"
"    barrier(CLK_GLOBAL_MEM_FENCE); \n"
"    \n"
"    temp = vector[row]; \n"
"    //eliminate column with index 'row' in parallel: \n"
"    for  (int elim = get_global_id(0); elim < row; elim += get_global_size(0)) \n"
"      vector[elim] -= temp * matrix[elim + row  * matrix_internal_rows]; \n"
"  } \n"
"   \n"
"}\n"
"__kernel void inplace_add(\n"
"          __global float * vec1,\n"
"          __global const float * vec2,\n"
"          unsigned int size) \n"
"{ \n"
"  for (unsigned int i = get_global_id(0); i < size; i += get_global_size(0))\n"
"    vec1[i] += vec2[i];\n"
"}\n"
"\n"
"__kernel void trans_unit_lower_triangular_substitute_inplace(\n"
"          __global const float * matrix,\n"
"          unsigned int matrix_rows,\n"
"          unsigned int matrix_cols,\n"
"          unsigned int matrix_internal_rows,\n"
"          unsigned int matrix_internal_cols,\n"
"          __global float * vector)\n"
"{\n"
"  float temp;\n"
"  for (int row = 0; row < matrix_rows; ++row)\n"
"  {\n"
"    barrier(CLK_GLOBAL_MEM_FENCE);\n"
"\n"
"    temp = vector[row];\n"
"\n"
"    for  (int elim = row + get_global_id(0) + 1; elim < matrix_rows; elim += get_global_size(0))\n"
"      vector[elim] -= temp * matrix[elim * matrix_internal_rows + row];\n"
"  }\n"
"}\n"
"\n"
"\n"
"__kernel void scaled_rank1_update(\n"
"          __global float * matrix,\n"
"          unsigned int matrix_rows,\n"
"          unsigned int matrix_cols,\n"
"          unsigned int matrix_internal_rows,\n"
"          unsigned int matrix_internal_cols,\n"
"          float val,\n"
"          __global const float * vector1,  \n"
"          __global const float * vector2) \n"
"{ \n"
"  float tmp;\n"
"  for (unsigned int row = get_global_id(0); row < matrix_rows; row += get_global_size(0))\n"
"  {\n"
"    tmp = val * vector1[row];\n"
"    for (unsigned int col = 0; col < matrix_cols; ++col)\n"
"      matrix[row + col*matrix_internal_rows] += tmp * vector2[col];\n"
"  }\n"
"}\n"
"__kernel void clear(\n"
"          __global float * vec,\n"
"          unsigned int size) \n"
"{ \n"
"  for (unsigned int i = get_global_id(0); i < size; i += get_global_size(0))\n"
"    vec[i] = 0;\n"
"}\n"
"__kernel void cpu_inplace_mult(\n"
"          __global float * vec,\n"
"          float factor, \n"
"          unsigned int size) \n"
"{ \n"
"  for (unsigned int i = get_global_id(0); i < size; i += get_global_size(0))\n"
"    vec[i] *= factor;\n"
"}\n"
"__kernel void inplace_mult(\n"
"          __global float * vec,\n"
"          __global const float * fac, \n"
"          unsigned int size) \n"
"{ \n"
"  float factor = *fac;\n"
"  for (unsigned int i = get_global_id(0); i < size; i += get_global_size(0))\n"
"    vec[i] *= factor;\n"
"}\n"
"__kernel void upper_triangular_substitute_inplace( \n"
"          __global const float * matrix, \n"
"          unsigned int matrix_rows,\n"
"          unsigned int matrix_cols,\n"
"          unsigned int matrix_internal_rows,\n"
"          unsigned int matrix_internal_cols,\n"
"          __global float * vector) \n"
"{ \n"
"  float temp; \n"
"  for (int row = matrix_rows-1; row > -1; --row) \n"
"  { \n"
"    barrier(CLK_GLOBAL_MEM_FENCE); \n"
"    if (get_global_id(0) == 0) \n"
"      vector[row] /= matrix[row + row*matrix_internal_rows]; \n"
" \n"
"    barrier(CLK_GLOBAL_MEM_FENCE); \n"
"    temp = vector[row]; \n"
"    //eliminate column with index 'row' in parallel: \n"
"    for  (int elim = get_global_id(0); elim < row; elim += get_global_size(0)) \n"
"      vector[elim] -= temp * matrix[elim + row  * matrix_internal_rows]; \n"
"  } \n"
"   \n"
"}\n"
;

const char * matrix_kernels_2 = 
"// file automatically generated - do not edit!\n"
"// matrix-matrix multiplication C = A^T * B^T\n"
"// matrix layouts: C...col_major, A...row_major, B...col_major\n"
"__kernel void prod_TT(\n"
"          __global const float * A,\n"
"          unsigned int A_rows,\n"
"          unsigned int A_cols,\n"
"          unsigned int A_internal_rows,\n"
"          unsigned int A_internal_cols,\n"
"          __global const float * B,  \n"
"          unsigned int B_rows,\n"
"          unsigned int B_cols,\n"
"          unsigned int B_internal_rows,\n"
"          unsigned int B_internal_cols,\n"
"          __global float * C,\n"
"          unsigned int C_rows,\n"
"          unsigned int C_cols,\n"
"          unsigned int C_internal_rows,\n"
"          unsigned int C_internal_cols,\n"
"          __local float * bufA,\n"
"          __local float * bufB) \n"
"{ \n"
"  int block_size = get_local_size(0);\n"
"  int row_block_id = get_group_id(0);\n"
"  int col_block_id = get_group_id(1);\n"
"  int row_thread_id = get_local_id(0);\n"
"  int col_thread_id = get_local_id(1);\n"
"  int aBegin = row_block_id * block_size;\n"
"  int aStep = block_size * A_internal_cols;\n"
"  int bBegin = col_block_id * block_size;\n"
"  int bStep = block_size * B_internal_rows;\n"
"  int block_num = A_rows / block_size;\n"
"  if (block_num * block_size != A_rows)\n"
"    ++block_num;\n"
"  float Csub = 0;\n"
"  int aOffset = row_thread_id + col_thread_id * A_internal_cols;\n"
"  int bOffset = row_thread_id * B_internal_rows + col_thread_id;\n"
"  for (int block = 0;\n"
"           block < block_num;\n"
"           ++block)\n"
"  {\n"
"    if (block * block_size + col_thread_id < A_rows && get_global_id(0) < A_cols)\n"
"      bufA[row_thread_id * block_size + col_thread_id] = A[aBegin + aOffset];\n"
"    else\n"
"      bufA[row_thread_id * block_size + col_thread_id] = 0;\n"
"    if ( (block * block_size + row_thread_id < B_cols) && get_global_id(1) < B_rows )\n"
"      bufB[row_thread_id * block_size + col_thread_id] = B[bBegin + bOffset]; \n"
"    else\n"
"      bufB[row_thread_id * block_size + col_thread_id] = 0;\n"
"    barrier(CLK_LOCAL_MEM_FENCE);\n"
"    for (int k = 0; k < block_size; ++k)\n"
"      Csub += bufA[row_thread_id * block_size + k] * bufB[k * block_size + col_thread_id];\n"
"    barrier(CLK_LOCAL_MEM_FENCE);\n"
"    aBegin += aStep;\n"
"    bBegin += bStep;\n"
"  }\n"
"  if (get_global_id(0) < A_cols && get_global_id(1) < B_rows)\n"
"    C[get_global_id(0) + get_global_id(1) * C_internal_rows] = Csub;\n"
"}\n"
"// file automatically generated - do not edit!\n"
"// matrix-matrix multiplication C = A^T * B\n"
"// matrix layouts: C...col_major, A...row_major, B...col_major\n"
"__kernel void prod_TA(\n"
"          __global const float * A,\n"
"          unsigned int A_rows,\n"
"          unsigned int A_cols,\n"
"          unsigned int A_internal_rows,\n"
"          unsigned int A_internal_cols,\n"
"          __global const float * B,  \n"
"          unsigned int B_rows,\n"
"          unsigned int B_cols,\n"
"          unsigned int B_internal_rows,\n"
"          unsigned int B_internal_cols,\n"
"          __global float * C,\n"
"          unsigned int C_rows,\n"
"          unsigned int C_cols,\n"
"          unsigned int C_internal_rows,\n"
"          unsigned int C_internal_cols,\n"
"          __local float * bufA,\n"
"          __local float * bufB) \n"
"{ \n"
"  int block_size = get_local_size(0);\n"
"  int row_block_id = get_group_id(0);\n"
"  int col_block_id = get_group_id(1);\n"
"  int row_thread_id = get_local_id(0);\n"
"  int col_thread_id = get_local_id(1);\n"
"  int aBegin = row_block_id * block_size;\n"
"  int aStep = block_size * A_internal_cols;\n"
"  int bBegin = col_block_id * block_size * B_internal_rows;\n"
"  int bStep = block_size;\n"
"  int block_num = A_rows / block_size;\n"
"  if (block_num * block_size != A_rows)\n"
"    ++block_num;\n"
"  float Csub = 0;\n"
"  int aOffset = row_thread_id + col_thread_id * A_internal_cols;\n"
"  int bOffset = row_thread_id + col_thread_id * B_internal_rows;\n"
"  for (int block = 0;\n"
"           block < block_num;\n"
"           ++block)\n"
"  {\n"
"    if (block * block_size + col_thread_id < A_rows && get_global_id(0) < A_cols)\n"
"      bufA[row_thread_id * block_size + col_thread_id] = A[aBegin + aOffset];\n"
"    else\n"
"      bufA[row_thread_id * block_size + col_thread_id] = 0;\n"
"    if ( (block * block_size + row_thread_id < B_rows) && get_global_id(1) < B_cols )\n"
"      bufB[row_thread_id * block_size + col_thread_id] = B[bBegin + bOffset]; \n"
"    else\n"
"      bufB[row_thread_id * block_size + col_thread_id] = 0;\n"
"    barrier(CLK_LOCAL_MEM_FENCE);\n"
"    for (int k = 0; k < block_size; ++k)\n"
"      Csub += bufA[row_thread_id * block_size + k] * bufB[k * block_size + col_thread_id];\n"
"    barrier(CLK_LOCAL_MEM_FENCE);\n"
"    aBegin += aStep;\n"
"    bBegin += bStep;\n"
"  }\n"
"  if (get_global_id(0) < A_cols && get_global_id(1) < B_cols)\n"
"    C[get_global_id(0) + get_global_id(1) * C_internal_rows] = Csub;\n"
"}\n"
"// file automatically generated - do not edit!\n"
"// matrix-matrix multiplication C = A * B\n"
"// matrix layouts: C...col_major, A...row_major, B...col_major\n"
"__kernel void prod_AA(\n"
"          __global const float * A,\n"
"          unsigned int A_rows,\n"
"          unsigned int A_cols,\n"
"          unsigned int A_internal_rows,\n"
"          unsigned int A_internal_cols,\n"
"          __global const float * B,  \n"
"          unsigned int B_rows,\n"
"          unsigned int B_cols,\n"
"          unsigned int B_internal_rows,\n"
"          unsigned int B_internal_cols,\n"
"          __global float * C,\n"
"          unsigned int C_rows,\n"
"          unsigned int C_cols,\n"
"          unsigned int C_internal_rows,\n"
"          unsigned int C_internal_cols,\n"
"          __local float * bufA,\n"
"          __local float * bufB) \n"
"{ \n"
"  int block_size = get_local_size(0);\n"
"  int row_block_id = get_group_id(0);\n"
"  int col_block_id = get_group_id(1);\n"
"  int row_thread_id = get_local_id(0);\n"
"  int col_thread_id = get_local_id(1);\n"
"  int aBegin = row_block_id * block_size * A_internal_cols;\n"
"  int aStep = block_size;\n"
"  int bBegin = col_block_id * block_size * B_internal_rows;\n"
"  int bStep = block_size;\n"
"  int block_num = A_cols / block_size;\n"
"  if (block_num * block_size != A_cols)\n"
"    ++block_num;\n"
"  float Csub = 0;\n"
"  int aOffset = row_thread_id * A_internal_cols + col_thread_id;\n"
"  int bOffset = row_thread_id + col_thread_id * B_internal_rows;\n"
"  for (int block = 0;\n"
"           block < block_num;\n"
"           ++block)\n"
"  {\n"
"    if (block * block_size + col_thread_id < A_cols && get_global_id(0) < A_rows)\n"
"      bufA[row_thread_id * block_size + col_thread_id] = A[aBegin + aOffset];\n"
"    else\n"
"      bufA[row_thread_id * block_size + col_thread_id] = 0;\n"
"    if ( (block * block_size + row_thread_id < B_rows) && get_global_id(1) < B_cols )\n"
"      bufB[row_thread_id * block_size + col_thread_id] = B[bBegin + bOffset]; \n"
"    else\n"
"      bufB[row_thread_id * block_size + col_thread_id] = 0;\n"
"    barrier(CLK_LOCAL_MEM_FENCE);\n"
"    for (int k = 0; k < block_size; ++k)\n"
"      Csub += bufA[row_thread_id * block_size + k] * bufB[k * block_size + col_thread_id];\n"
"    barrier(CLK_LOCAL_MEM_FENCE);\n"
"    aBegin += aStep;\n"
"    bBegin += bStep;\n"
"  }\n"
"  if (get_global_id(0) < A_rows && get_global_id(1) < B_cols)\n"
"    C[get_global_id(0) + get_global_id(1) * C_internal_rows] = Csub;\n"
"}\n"
"// file automatically generated - do not edit!\n"
"// matrix-matrix multiplication C = A * B^T\n"
"// matrix layouts: C...col_major, A...row_major, B...col_major\n"
"__kernel void prod_AT(\n"
"          __global const float * A,\n"
"          unsigned int A_rows,\n"
"          unsigned int A_cols,\n"
"          unsigned int A_internal_rows,\n"
"          unsigned int A_internal_cols,\n"
"          __global const float * B,  \n"
"          unsigned int B_rows,\n"
"          unsigned int B_cols,\n"
"          unsigned int B_internal_rows,\n"
"          unsigned int B_internal_cols,\n"
"          __global float * C,\n"
"          unsigned int C_rows,\n"
"          unsigned int C_cols,\n"
"          unsigned int C_internal_rows,\n"
"          unsigned int C_internal_cols,\n"
"          __local float * bufA,\n"
"          __local float * bufB) \n"
"{ \n"
"  int block_size = get_local_size(0);\n"
"  int row_block_id = get_group_id(0);\n"
"  int col_block_id = get_group_id(1);\n"
"  int row_thread_id = get_local_id(0);\n"
"  int col_thread_id = get_local_id(1);\n"
"  int aBegin = row_block_id * block_size * A_internal_cols;\n"
"  int aStep = block_size;\n"
"  int bBegin = col_block_id * block_size;\n"
"  int bStep = block_size * B_internal_rows;\n"
"  int block_num = A_cols / block_size;\n"
"  if (block_num * block_size != A_cols)\n"
"    ++block_num;\n"
"  float Csub = 0;\n"
"  int aOffset = row_thread_id * A_internal_cols + col_thread_id;\n"
"  int bOffset = row_thread_id * B_internal_rows + col_thread_id;\n"
"  for (int block = 0;\n"
"           block < block_num;\n"
"           ++block)\n"
"  {\n"
"    if (block * block_size + col_thread_id < A_cols && get_global_id(0) < A_rows)\n"
"      bufA[row_thread_id * block_size + col_thread_id] = A[aBegin + aOffset];\n"
"    else\n"
"      bufA[row_thread_id * block_size + col_thread_id] = 0;\n"
"    if ( (block * block_size + row_thread_id < B_cols) && get_global_id(1) < B_rows )\n"
"      bufB[row_thread_id * block_size + col_thread_id] = B[bBegin + bOffset]; \n"
"    else\n"
"      bufB[row_thread_id * block_size + col_thread_id] = 0;\n"
"    barrier(CLK_LOCAL_MEM_FENCE);\n"
"    for (int k = 0; k < block_size; ++k)\n"
"      Csub += bufA[row_thread_id * block_size + k] * bufB[k * block_size + col_thread_id];\n"
"    barrier(CLK_LOCAL_MEM_FENCE);\n"
"    aBegin += aStep;\n"
"    bBegin += bStep;\n"
"  }\n"
"  if (get_global_id(0) < A_rows && get_global_id(1) < B_rows)\n"
"    C[get_global_id(0) + get_global_id(1) * C_internal_rows] = Csub;\n"
"}\n"
"__kernel void unit_lower_triangular_substitute_inplace(\n"
"          __global const float * matrix,\n"
"          unsigned int matrix_rows,\n"
"          unsigned int matrix_cols,\n"
"          unsigned int matrix_internal_rows,\n"
"          unsigned int matrix_internal_cols,\n"
"          __global float * vector)\n"
"{\n"
"  float temp;\n"
"  for (int row = 0; row < matrix_rows; ++row)\n"
"  {\n"
"    barrier(CLK_GLOBAL_MEM_FENCE);\n"
"    temp = vector[row];\n"
"    for  (int elim = row + get_global_id(0) + 1; elim < matrix_rows; elim += get_global_size(0))\n"
"      vector[elim] -= temp * matrix[row * matrix_internal_rows + elim];\n"
"  }\n"
"}\n"
"__kernel void inplace_sub(\n"
"          __global float * vec1,\n"
"          __global const float * vec2,\n"
"          unsigned int size) \n"
"{ \n"
"  for (unsigned int i = get_global_id(0); i < size; i += get_global_size(0))\n"
"    vec1[i] -= vec2[i];\n"
"}\n"
"__kernel void lower_triangular_substitute_inplace(\n"
"          __global const float * matrix,\n"
"          unsigned int matrix_rows,\n"
"          unsigned int matrix_cols,\n"
"          unsigned int matrix_internal_rows,\n"
"          unsigned int matrix_internal_cols,\n"
"          __global float * vector)\n"
"{\n"
"  float temp;\n"
"  for (int row = 0; row < matrix_rows; ++row)\n"
"  {\n"
"    barrier(CLK_GLOBAL_MEM_FENCE);\n"
"    if (get_global_id(0) == 0)\n"
"      vector[row] /= matrix[row+row*matrix_internal_rows];\n"
"    barrier(CLK_GLOBAL_MEM_FENCE);\n"
"    temp = vector[row];\n"
"    for  (int elim = row + get_global_id(0) + 1; elim < matrix_rows; elim += get_global_size(0))\n"
"      vector[elim] -= temp * matrix[row * matrix_internal_rows + elim];\n"
"  }\n"
"}\n"
"__kernel void trans_vec_mul(\n"
"          __global const float * matrix,\n"
"          unsigned int matrix_rows,\n"
"          unsigned int matrix_cols,\n"
"          unsigned int matrix_internal_rows,\n"
"          unsigned int matrix_internal_cols,\n"
"          __global const float * vector,  \n"
"          __global float * result) \n"
"{ \n"
"  //row and col indicate indices within transposed matrix\n"
"  for (unsigned int row = get_global_id(0); row < matrix_cols; row += get_global_size(0))\n"
"  {\n"
"    float dot_prod2 = 0.0f;\n"
"    for (unsigned int col = 0; col < matrix_rows; ++col)\n"
"      dot_prod2 += matrix[row * matrix_internal_rows + col] * vector[col];\n"
"    result[row] = dot_prod2;\n"
"  }\n"
"}\n"
"//perform a rank-1 update of the matrix, i.e. A += x * x^T\n"
"__kernel void rank1_update(\n"
"          __global float * matrix,\n"
"          unsigned int matrix_rows,\n"
"          unsigned int matrix_cols,\n"
"          unsigned int matrix_internal_rows,\n"
"          unsigned int matrix_internal_cols,\n"
"          __global const float * vector1,  \n"
"          __global const float * vector2) \n"
"{ \n"
"  float tmp;\n"
"  for (unsigned int row= get_global_id(0); row < matrix_rows; row += get_global_size(0))\n"
"  {\n"
"    tmp = vector1[row];\n"
"    for (unsigned int col = 0; col < matrix_cols; ++col)\n"
"      matrix[row + col * matrix_internal_rows] += tmp * vector2[col];\n"
"  }\n"
"}\n";
